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We study the ground-state properties of the double-chain Hubbard model coupled with ferro- 
magnetic exchange interaction by using the weak-coupling theory, density-matrix renormalization 
group technique, and Lanczos exact-diagonalization method. We determine the ground-state phase 
diagram in the parameter space of the ferromagnetic exchange interaction and band filling. We find 
that, in high electron density regime, the spin gap opens and the spin-singlet d^-wave-like pairing 
correlation is most dominant, whereas in low electron density regime, the fully-polarized ferromag- 
netic state is stabilized where the spin-triplet p^-wave-like pairing correlation is most dominant. 

PACS numbers: 71.10.Fd, 71.10.Hf, 74.20.Mn, 74.20.Rp 



I. INTRODUCTION 



Ferromagnetism in itinerant electron systems has in- 
creasingly been understood since the Hubbard model 
was introduced in 1963.^2 It is known that subtle in- 
terplay between quantum many-body effects and spin- 
independent Coulomb interactions plays a crucial role in 
generating ferromagnetic orders in some solidsj 3 - a vari- 
ety of origins such as Nagaoka-Thouless mechanism, 4,5 
flat-band^ orbital degenerac y 8 ' 9 ' 10 three-site ring ex- 
change interaction ) 11 ' 12 etc., have so far been proposed. 
In addition to ferromagnetism, (possibly) spin-triplet su- 
perconductivity was recently discovered in the metallic 
ferromagnets UGe 2f ^ URhGe^ and ZrZna-i 5 - Conse- 
quently, the relation between superconductivity and fer- 
romagnetism has been of special interest in the field of 
strongly correlated systems. From the theoretical point 
of view, the occurrence of superconductivity in ferromag- 
netic materials is naturally explained by the formation 
of Cooper pairs with parallel spins, namely, spin-triplet 
pairs J£ 

Among the origins of ferromagnetism mentioned 
above, only the three-site ring exchange interaction acts 
on a couple of electrons. It yields a ferromagnetic spin 
correlation, which in turn produces an attractive effect 
between them. One may easily imagine that a spin-triplet 
superconductivity is realized if the attractive interaction 
between electrons can survive against the other effects. 
Recently, we have confirmed that this mechanism actu- 
ally works in a fairly simple correlated electron system; 
the system consists of two Hubbard chains coupled with 
zigzag bonds and has a unique structure of hopping intc- 
grals i 17 ' 18 ' 19 In this model, the spin-triplet pairing of elec- 
trons occurs between the inter-chain neighboring sites. If 
the ferromagnetic correlation between the two chains are 



essential for the spin-triplet superconductivity, we may 
be allowed to mimic our original model by a double-chain 
Hubbard model coupled with ferromagnetic exchange in- 
teraction. This new model is much easier to analyze than 
the original one due to the reduction of quantum fluctu- 
ations. Therefore, the introduction of this model will en- 
able us to investigate the spin-triplet superconductivity 
in more detail for wide range of parameters. 

In this paper, we thus study the double-chain Hubbard 
model coupled with ferromagnetic exchange interaction. 
We use the weak-coupling bosonization/renormalization 
group (RG) analyses^ - the density-matrix renormaliza- 
tion group (DMRG) technique, and the Lanczos exact- 
diagonalization method^ We thereby determine the 
ground-state phase diagram: we find that, in the high 
electron density regions, the spin gap opens and the 
d^y-wave-like pairing occurs, while in the low electron 
density regions, the system is in the metallic state with 
full spin polarization, where the pj,-wave-like spin-triplet 
pairing correlation becomes the most dominant. We note 
that this model can be regarded as a single-chain model 
with two degenerate orbitals if we assume the ferromag- 
netic exchange interaction to be identified with the intra- 
atomic Hund's rule coupling. This single-chain model has 
so far been studied both analyticall y 23 ' 24 ' 25 and numer- 
icall y 26 ' 27 . It has been proposed that the Haldane gap 
state is realized at half filling. It has also been pointed 
out that the system remains gapful for low hole doping 
regions i^ 8 . 

Our paper is organized as follows. In Sec. |Hl we de- 
fine the double-chain Hubbard model coupled with fer- 
romagnetic exchange interaction. In Sec. IIIH we analyze 
the model using the weak-coupling theory and derive the 
pairing order parameter. In Sec. IIV1 we calculate several 
quantities with the DMRG and Lanczos methods and 
present the calculated results. We also check if the results 
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FIG. 1: Schematic representation of the double-chain Hub- 
bard model coupled with ferromagnetic exchange interaction. 
No hopping process is allowed between the two chains and no 
exchange interactions exist along the chains. 



are consistent with the weak-coupling results. Section fVl 
contains summary and conclusions. 



II. MODEL 

We study the double-chain Hubbard model coupled 
with ferromagnetic exchange interaction (see Fig. [1]), the 
Hamiltonian of which is defined by 

H= 1 12 {4 x ,r y , a Cra: + l,r y ,a + H.C.J 

+U ^2 1-^,^,4. 

~j12 Sr -' 1 ■ Sr -> 2 w 



where cj^ r a {c r!C ,ry,a) is the creation (annihilation) op- 
erator of an electron with spin a (=T>4) at site r x on 
leg r y (= 1,2), n rxfrv>a = 4„r v , < 7 C n»>r v ,ff is the density 
operator, and S ra ., r is the spin-i operator, t is the 
nearest-neighbor hopping integral along the chain, U is 
the onsite Coulomb interaction, and J(> 0) is the ferro- 
magnetic exchange interaction between two sites on each 
rung. Note that no hopping process is allowed between 
the two chains and thus the two bands are degenerate. 



III. WEAK-COUPLING THEORY 

We first consider the weak-coupling regime where only 
the low-energy excitations near the Fermi points are cru- 
cial. Thus far, the weak-coupling theory has been well- 
constructed for two-band models i 2 ^^^^^^^^^ 
We further develop the approach to analyze the Hamil- 
tonian ([1]). Assuming a linearization of the dispersion 
relations in the vicinity of the Fermi level, we introduce 
the field operators of right- and left-going electrons as 
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where Cp a is the Fourier transform of combined operator 

(Cr x ,\,<r ± Cr x ,l,<r)/y/% for the right-going (p = +) and 
left-going (p = — ) electrons. L is the chain length where 
the lattice spacing is set to be unity. Using these field 



operators, the Hamiltonian can be rewritten as a sum of 
the linearized kinetic energy and interaction terms. We 
thus obtain 



H 



dxHn 



dxHi 



H = v F Y, ^la,C ( -^fa ) 



Hl = 1 1212 'all^la^-p,-^^-^- 



P,<r C 



4 12 12 '92l^i,a^- P ,-a,c, 2 ^-P,-<y,C^P,<yX 3 

ll212 'dl'^'l-Xi^P^^C^-P^Xs- (3) 
v,° c 

where e = £1(3 an d e = ClC2- The primed summa- 
tion over Q(i = 1,2,3,4) is restricted by a relation 
C1C2C3C4 = +1; which comes from the momentum con- 
servation condition in the transverse components. The 
coupling constants and g^x are related to the origi- 
nal parameters in the Hamiltonian (JTJ) : 
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For the SU(2) symmetric case, we can choose 
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(6) 
(7) 
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In this Section III, we consider the case away from half 
filling. Hence, the Umklapp term g^ gives no contribu- 
tion and the Fermi velocity renormalization due to the 
forward-scattering term 54 may be neglected. 



A. Bosonization 

Using the Abelian bosonization method^ we intro- 
duce eight chiral bosonic fields <$P where \i refers to the 
charge (p) and spin (a) sectors; meanwhile, r refers to 
the even (+) and odd (— ) sectors. The bosonic fields 
satisfy the commutation relations [<^ r (x), 4>% r '{ x ')] = 
±i(7r/4)sgn(a; - x')S^^S r y and [</>£ r (x),</>~ /r ,(aO] = 
j(7r/4)5 Mj/J '(5 r . r '. We then define a new set of chiral 
bosonic fields as 



(9) 



where p = ±, s — ±, and £ = ±. The chiral bosons obey 
the commutation relations \phip tSl £(x),phi PtS > l £'(x')] = 
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±i(w/4)sga{x-x')Sp t p>S Sl gi and [(p +<s ^(x), 0_ )S /,^(a;')] = 
i(w/4:)S StS tS^^/. The field operators of the right-moving 
and left-moving electrons ([2]) are then written as 



\j2ma 



exp (ipkpx + ip4> PiS ^) 



(10) 



where s = + (— ) for a —\ (J,). The Majorana fermions 
rj a ^, known as Klein factors, are introduced to ensure 
the proper anticommutation relations between fermion 
fields with different band and spin indices. They obey 
{rja-Xi Va-'X'} ~ ^T,a-'^CX'- ^ ^ s generally more convenient 
to trade the chiral boson fields pairwise for a conventional 
bosonic phase field (j> and its dual field 9, so that we also 
introduce the bosonic fields given by 
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' fj,,r " 



rfl,r 



(11) 

(12) 



where the operator H tir (x) = OxO^/ir is a canonical con- 
jugate variable to 4>^ r , which satisfies [4>^ r {x), II Mr (x')] = 
iS(x - x')8 tl ^5 r y. 

Now, we can rewrite the noninteracting term of the 
Hamiltonian as 



V 



2tt ^ 

fi—p,(T r—± 



^2 Vp^^rf + {9x4>n,r) 



(13) 



On the other hand, the interacting term is more com- 
plicated, containing many products of the Klein fac- 
tors such as f = .+ r ll-+ r lT --^T h a = Vcr.+Va-.-, an d 
Iiq = rjf^rji^. However, it is known that their eigenval- 
ues are T = ±1, h a = ±i, and h'^ = ±i.~ Thus, if we 
adopt the following convention T = +1, h a = i, h's = i(, 
the interacting term is reduced to 



H, = 



E f^^A^r 



E 



2^2 [{9i± - 9 2 ± ) cos2(?!) p _ cos2( ? (> ._ 



+ (ffix - 9 2 ±) cos 26 p 
+gf± cos 2tj) a+ cos 2<f> a 

~9\X cos 2 ®p- cos 2(/) cr - 



cos 29 „- 

+ g^ cos 20 p _ cos 2(j) a+ 
+ g±£ cos 2<f) a+ cos 29 „- 
g^ cos 29 p _ cos 20 ff _ + cos 20 p _ cos 2# CT _] 

(14) 



with 
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9 P - = 


~9tt 


+■ 2.9 2 \ + - 


-9iT - 
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(16) 


9*+ = 


-9tl 


- 2ffrT 






(17) 


9a- = 


~9tt 


+ 2 5l T- 






(18) 



If we use the notation with g^^r, Ho can also be written 



^0 - E if K,r (5,^,,) 2 + (^0^ r ) 2 l (19) 



with the critical exponents 



27TVF - ffp.r 



27TW F + g^, r 

and the renormalized Fermi velocity 



(20) 



(21) 



B. Renormalization-group analysis 

By treating the interaction perturbatively, we derive 
the renormalization group (RG) equations in the one- 
loop level as follows: 



dg+ 



1_L 



OS 



= -2te + ) -2g-+g-+ 



-2 (at) +2.9 1 X.9 2 



2X 



(22) 



as 



-2 ffl +.9++ - 2<? 2 +g++ - 4 9l +g lx 

-2 fc + .9 2 + r + .9iT5 2 Y) (23) 



dS 



-2(g- 1 +) I -2(g^y 
-2g-+g-+ - 2g+£g^ 



(24) 



9gt± 
dS 



-2 (2g+- - g--) g++ + 2g+- g++ 
' 2 9i±92± - 2 9t±9t± (25) 



dg+ 



2 _L 



dS 



- fe + ) 2 - fc + ) 2 - fe + ) 2 + feT) 2 (26) 



= - 2 9il9il 2 9^92± + 2 9 2 Utl (27) 



' hjl - - 2 9i±9tZ + 2 9tt9iZ - 2 92±9tl (28) 



dS 



d 92± 

dS 



~(9^) -(9tL) 
+ (92±) 2 - (92±) 2 



(29) 



where S = In I /(ivvf) is the RG time and I is the scaling 
quantity. We note that the couplings g^± and g%x play 
important roles here. In general, those couplings are ir- 
relevant when the hopping processes between two chains 
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are relevan t 24 i 32 ' 33 ' 34 i 36 . However, there is no transverse 
hopping term in our model, so that we have to start from 
the picture of degenerate bands. 

By solving the RG equations (|2"2")) - (l2l)l) , we obtain the 
relations 



9n >32_l ' 9i± 7 

92± 



-ati. » Sax » 9tt » 



(30) 



and g]^ = g 2 , + = for U > J. In this case, the pa- 



rameters Kp-i K a+ , and K a - are scaled as K p 



0, 



K, 



0, and K a - — > oo. We therefore can simplify the 



second term of the phase Hamiltonian (fl4|) as 



(f ix - 9%i ) cos 2 0p- cos 2(j) a 



+ffix cos 20<t+ cos 2(b a 
+9ii cos 2(f) IJ+ cos 20 CT 



- g^j_ cos 2</> p _ cos 20 CT + 
SzjT cos 2</>p- cos 20 CT _ . 



(31) 



Note that the phase variables 9 P - and 9„ + are not in- 
cluded in this term. Taking g^ > into account, wc 
find that the fields (b p - an d (f> a + are locked at (<f> p -) = 



%I\ + 7r/ 2 and 



= 5(Ji + 1) + 7r/3 where /„ are 



integers. 

Let us then turn to the a— mode that remains to be 
studied. The effective Hamiltonian of the a— mode takes 
the form 



= dx 



+g<p cos 20cr_ + cos 29a-] 



(32) 



where g^ = g+ ± - g 2± - and g g = g 2± - g 1± . We 
here adopt a set of the variables (<j) P - , 4> a + ) = {^h > "7 2 + 
71/2); however, it leads to a physically equivalent result if 
we chose another set ((b p _,(b a+ ) = (n/2 + nli,irl 2 )- For 
the Hamiltonian (|32|) . there are three RG equations as 
follows: 



dl 
dy4, 

dl 
dye 

dl 



Ve - K-Vl 
{2-2K tT -)y <j) 

(2-2K-±)y e . 



(33) 
(34) 
(35) 



where y a = g a /(-KVF) with a = 6,<b. Since \ye\ 3> \y^>\ 



and K„- > 1, the parameters are renormalized to ye 
°°! y<t> ~ * 0, and Kcr- — > 00. As a results, we find that 

f/i+7r/ 2 ,^+ = f(/i + 



three modes arc locked as . 



1) + 71-/3, and 9 a - = + 1) + 7r/ 4 . Therefore, there is 
a gap in all the spin excitations. 



C. Order parameters 

In order to determine the dominant correlation, we in- 
troduce the order parameters of the superconducting cor- 
relations for four kinds of pairing symmetry, which are 
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FIG. 2: Symmetry of the pairing state in the double-chain 
Hubbard model. The sign of the order parameter is also 
shown. 



shown in Fig. [2j They are given in terms of the phase 
fields as follows: 

Osi = ^hC^-Pd-C 
p=±.C=± 
2 

~ — e l9p+ [cos 6 a - cos 4>,j + cos 6 P - 
■kcl 

—i sin 9c- sin (b a + sin (b p -] 

On = C^WP-pd -i 

p=±.C=± 



(36) 



S 2 = 



— e l p+ [i sin 6* CT _ cos cb a+ cos </>„_ 
na 

— cos 9 a - sin (b a+ sin <b p - ] 

Y CM>,u^-p,i,-c 

p=±,C=± 
2i 

— e l6,p+ \i sin 9 a - sin 6 a + cos <b - 



(37) 



Ot2 = 



■ COS ti a - COS ( 



sm< 



(38) 



p=±,C=± 

2i ;a , „ . , 



cos tier- sm< 



-i sm t' (T _ cos CT + sin < 



cos O - 



(39) 



One can easily find that the inter-chain diagonal singlet 
pairing state, which corresponds to the order parameter 
Os2, is only of the quasi-ordering superconducting in- 
stability when (b p - and 9 a - are locked. The asymptotic 
behavior of the S2 correlation function is given by 



C4(r)0 s2 (0) 



1 



A/2K P 



(40) 



and the other interchain pairing correlations decay expo- 
nentially. We also mention that all the intrachain pair- 
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ing correlations decay exponentially due to the locked 
modes. If we assume a S2 ground state 0' S2 |0) and ap- 
ply the interacting term of Hamiltonian (fT4"]) to the state, 
we obtain 

HjO\ 2 10) = fo-f - g+Z - ff 2 T + .g 2 + D0^ 2 |0) . (41) 

Since g^ < 0, g±£ > 0, g^jj > 0, and g^T < 0, we find 
that the S2 pairing state can gain much energy. 

The other possible order is the 4fcp charge-density wave 
(CDW) correlation, of which the order parameter is given 
by 

4 fe F = ^P,<T,+^t,<7',-^-P^'-^-P,^ + 

~ -5—^ cos (4fcFX + 20 p+ ) (cos + cos 2d> a - ) , 

(42) 

where the cos 2<j} a+ term comes from the cases where the 
spins are parallel a = a' , and the cos2</v_ term comes 
from the cases where the spins are antiparallcl a = — a' . 
The latter component cos 2(j) a ^ decays exponentially be- 
cause the field is locked. Thus, the 4fcp-CDW corre- 
lation shows a power-low behavior like 



T| 1 1 I I I I 1 1| 1 1 — I I I I ll| 1 1 I I I I 1 1| 

hx — x x — Mmonin«»o-nn * — ^ 



((4 F (r)0 4feF ) 



1 



(43) 



By comparing Eq.(|35]) with Eq.fgjj]), we find that the 
d^-wave-likc S2 pairing correlation is most dominant for 
K p+ > 0.5, whereas the 4/cf-CDW is most dominant for 
K p+ < 0.5. This is the same as that of the standard 
two-band model. 33 



IV. NUMERICAL RESULTS 

Next, we turn to the intermediate to strong coupling 
regime. We employ the Lanczos and DMRG methods 
to obtain energies and physical quantities in the ground 
state and low-lying excited states. In order to carry out 
our calculations, we consider N (= N^+N^) electrons in a 
system with length L (containing 2L sites) . The electron 
density is given by n — N/L. Note that the number of 
electrons must be taken as N — 4/ with Z (> 1) being 
an integer to maintain the total spin of the unpolarizcd 
ground state as S = 0. 

For static quantities, we use the DMRG method with 
applying the open-end boundary conditions for precise 
calculations. We study systems with length up to L = 
128 and keep up to m w 2400 density-matrix eigen- 
states in the DMRG procedure. In this way, the dis- 
carded weights are typically of the order 10~ 8 ~ 10~ 7 
and the ground-state energy is obtained in the accuracy 
of ~ 10 _7 i. All the calculated energies are extrapolated 
to the limit m — > 00. For dynamical quantities, we use 
the Lanczos method for small clusters with the periodic 
boundary conditions. The system size is assumed as 
L = 8, i.e., 8x2 ladder. 




■ ^' " tja-ij-Darj-^ — 

100 



FIG. 3: (Color online) Calculated values of the total-spin 
quantum number S as a function of J for U — 20 (circles), 40 
(squares), and 00 (crosses). L = 48 and n = 0.5 are assumed. 



A. Spin polarization 

Approaching from the strong-coupling regime U t, 
we may anticipate the presence of the fully-polarized fer- 
romagnetic state. When the two chains are uncoupled, 
i.e., J = 0, the state can be interpreted as the Nagaoka 
state. Generally, the appearance of the Nagaoka state 
is limited to large-[7 and low-ro range. However, the 
fully-polarized region would be spread into smaller-?/ and 
higher-n range if a finite J is introduced. This is because 
the polarized electrons can gain the kinetic energy with- 
out loss of the exchange interaction between the chains. 

Let us start with investigating how the ferromagnetic 
phase appears in the parameter space (U,J). We can 
find it by calculating the expectation value of total-spin 
operator S in the ground state, which is defined by 



<S 2 ) = ^(S l -S,) = ^ + l) 



(44) 



For a fully-polarized state, one will obtain the total-spin 
quantum number S — S max = N/2, i.e., S/S max = 1. In 
Fig. [3J we show the total spin S normalized with respect 
to Smax as a function of J for several values of U. The 
system size and filling are fixed at L = 48 and n = 0.5, 
respectively. We calculate S for systems with length L — 
24, 36, and 48, and confirm that the size dependence is 
negligible. 

At U = 20, with increasing J, we find a transition 
from the paramagnetic state to the fully-polarized state 
at J c i ~ 4.5. Moreover, we find that the system goes back 
to the paramagnetic state at J C 2 ~ 8.5. This is because 
the formation of the local spin-singlet bound states gains 
more energy for very large values of J since J increases 
the gain in the spin-singlet binding energy as shown in 
Sec lIV B"3l We also find that the region with the full spin 
polarization broadens with increasing U. In the limit of 
U = 00, J c i — * and J c i — * 00. We note that both of 
the transitions are discontinuous, i.e., of the first-order 
in the thermodynamic limit. Thus, the critical transition 
points can be determined in the parameter space (U, J), 
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near half filling 



FIG. 4: (Color online) A C (L) as a function of the inverse sys- 
tem length 1/L. Lines are the fits with the 2nd-order poly- 
nomials. 



which will be given as a ground-state phase diagram in 
Sec lIVFl 



B. Excitation gaps 

To ascertain the presence of the lowest excitation gap 
in the charge, spin, and pairing sectors, we calculate the 
charge gap, spin gap, and binding energy in the ther- 
modynamic limit. We study several chains with lengths 
up to L = 64 and then perform the finite-size scaling 
analysis based on the size-dependence of each quantity. 



1. Charge gap 
The charge gap A c is defined by 
A c = lim A C {L), (45) 

L— »oo 

A C (L) = £°(iV t + 1, N t + 1) + El (JV T -l,N x - I) 
-2E° L (N e ) 

where E^ (iVj , 2Vj ) denotes the ground-state energy of a 
chain of length L with N-\ spin-up and N± spin-down 
electrons. In Fig. 0] we plot the charge gap A C (L) as 
a function of the inverse system length 1/L for several 
parameter sets. 

At half filling (n = I), we can easily expect the system 
to be a Mott insulator for U > 0. However, we also find 
the system is insulating even for U — if J is finite. This 
is associated with the fact that a localized spin-triplet 
pair is formed on each rung. If an electron is add to the 
rung, an effective on-rung repulsive interaction U c g = 
U + i acts. Thus, we obtain the effective Hamiltonian 



u - 



J 



where c rxj(T = cv i i<T c rxi 2,<T is the annihilation operator 
of an electron with spin a at rung r x , n rx . = cj. ^c^^ 
is the number operator, S* = nr m ,l ~~ n r m ,i is the z- 
component of total spin, and = TJ T rj . with the spin- 
triplet operator T Ti = c Tat i < a c r*&,<' or ^( £ V«.i,o- c r-<. J 2,ff + 

C-r^ ,l,CTCr x ,2,cr ■ 

Let us then turn to the case away from half filling. 
There are two points to be noted here: one is whether the 
system is insulating at quarter filling (n = 0.5) as is the 
case of the Hubbard and t—J ladder; the other is whether 
the phase separation occurs in the fully-polarized state. 
As seen in Fig. 01 all the results for n = 0.5 are smoothly 
extrapolated to zero in the thermodynamic limit 1/L — » 
0. At U = oo and J = 4, the system is in the fully- 
polarized state. We therefore conclude that the system 
is metallic in the entire region except n = 1. 



2. Spin gap 



The spin gap A s is defined by 
A s = lim A S {L), 

L — >co 

A S (L) = E° L (Ni + 1,^-1)- EKN^NJ 



(47) 



In Fig. [5l we plot the spin gap A S (L) as a function of 
inverse system length 1/L at (a) n = 1 and (b) n = 0.75. 
For U = 4, n = 1 and U = 1, n = 0.75, values of A S (L) 
seem to be smoothly extrapolated to zero when 1/L—* 0, 
which is in contrast to the weak-coupling analysis. For 
the other cases, however, the extrapolated lines reach 
zero at finite values of 1/L. This is unphysical and so 
it may be a good guess that A S (L) is fairly flat around 
1/L = 0, as seen in the results for U — 4,n — 1 and 
U = l.n = 0.75. Unfortunately, we cannot treat the 
large enough systems to confirm if this is the case. Nev- 
ertheless, if we assume that the 1 / L-dependence of A s (L) 
reflects the spinon band structure near the Fermi point, 
the flat behavior of A S (L) may rather imply the presence 
of a spin gapful state. Since we cannot detect the spin gap 
less than ~ I0~ 6 within the present calculations, the ex- 
istence of a very small spin gap, i.e., A(L — > oo) < I0~ 6 , 
is conceivable. 

To determine if the spin gap is present or absent, we 
also calculate the equal-time spin-spin correlation func- 
tion 



S(\r x 



) ~ (^r T ,r„ ^r' .r' ) (^r T ,r„ ) .r'„ ) 1 



where r ^ is the z component of the spin operator of 
an electron at site r x on rung r y . In Fig. we present 



0.05 0.1 
1/L 



0.05 0.1 

1/L 




FIG. 5: (Color online) Upper panels: A„(L) as a function of 
inverse system size for (a) n = 1 and (b) n = 0.75. Solid lines 
are the polynomial fits to the data for the finite-size scaling 
analysis. Lower panels: Semilogarithmic plot of the magni- 
tude of the spin-spin correlation function S(]r x — r' x \, 1, 2') at 
J — 1 for (c) n = 1 and (d) n = 0.75. The data are fitted 

with a function S(\r x — r' x \, 1, 2) ~ exp 




n — i — i — 1 — i — 1 — r 



U=2Q 



(c) 



Fully-polarized \ 
region \ 



10 



FIG. 6: (Color online) Upper panel: A^(L) as a function 
of inverse system length 1/L at n = 0.5. Solid lines are 
the polynomial fits to the data for finite-size scaling analy- 
sis. Solid (empty) symbols denote the binding energy of elec- 
trons (holes). Lower panels: Values of Ag(I) extrapolated 
to 1/L — > plotted as a function of J for (b) U — 4 abd (c) 
U = 10. 



3. Bin 



energy 



a semilogarithmic plot of l^dr^ — r' x \, 1, 2)| as a function 
of the distance \r x — r' x \ for (c) n = 1 and (d) n = 0.75. 
Note that the long-range behavior of I^Qr^ — r' x \, 1, 1)| is 
almost the same as that of\S(\r x — r' x \,l,2)\. The results 
can be fitted with a function exp(— \r x — r' x \/^) and thus 
the exponential decay of the correlation functions is con- 
firmed for all the parameter sets. The correlation lengths 
are estimated as £ = 10.35(9.74) for U = 1(4) at n = 1; 
£ = 14.29(12.99) for U = 1(4) at n = 0.75. According to 
the Tomonaga-Luttinger liquid theory, the spin-spin cor- 
relation can decay exponentially only when there is a gap 
in all the magnetic excitations. Consequently, it would 
be rather appropriate to conclude that there exists quite 
small spin gap (< 10 -6 ). The correlation lengths seem to 
be much longer than those of other standard spin-gapped 
systems, e.g.,£ = 3.19 in the two-leg isotropic Heisenberg 
system. However, it has been found that in the zigzag 
Heisenberg chain, the correlation lengths increase rapidly 
with decreasing binding energy^ Thus, the very large 
values of £ may reflect an exponentially small spin gap. 



The binding energy A B is defined by 

A| = lim A±(L), (49) 

L — >oo 

A±(L) - El (7V T ± 1, N l ± 1) + E° L (N h Ni) 

where the + (— ) sign corresponds to the binding energy 
of electrons (holes). In Fig. [S] (a), we plot the binding 
energy A^ (L) as a function of the inverse system length 
1/L at n = 0.5 for several sets of parameters. For all the 
parameter sets, |A^(L)| is found to decrease monoton- 
ically as a function of 1/L, so that we can extrapolate 
|A^(L)| to the thermodynamic limit systematically. We 
perform a least-squares fit of |A^(L)| to a polynomial in 
1/L and obtain the extrapolated values. We note that 
the binding energies of electrons and holes are extrapo- 
lated to the same value at 1/L — > 0. 

In Fig.[S](b), the extrapolated values of A^ for U = 4 
as a function of J are shown. When U is small, the sys- 
tem is expected to be in the spin-singlet superconducting 
phase, as discussed in Sec. IIIII Hence, the binding en- 
ergy is determined by an energy of the spin-singlet bound 
state. In analogy with the Haldane gap, we expect a scal- 
ing A^ cx Jt/U 2 from perturbation assuming the double 
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FIG. 7: (Color online) Log-log plot of the pair-correlation 
functions D a (\r x — r x \) calculated for L — 128 at (a) n = 
0.825, U = 10, J = 2, (b) n = 0.5, U = 10, J = 4, and (c) 
n = 0.25, U = 10, J = 5. 



occupancy is sufficiently excluded. The origin of this en- 
ergy gain is also interpreted as the double exchange inter- 
action in our model. For small J, this estimation seems 
not to work well. However, the double occupancy is in- 
creasingly excluded with increasing J and thus |A^(L)| 
increases linearly with J at J > 3. fn Fig. [5] (c), the ex- 
trapolated values of for U = 20 as a function of J are 
shown. Since the double occupancy is almost excluded, 
|Ap| increases linearly with J even if J is small. It is 
notable that the binding energy is strongly suppressed in 
the fully-polarized phase. This phase is regarded as the 
spin-triplet superconducting one, where the binding en- 
ergy is determined by an energy of the spin-triplet bound 
state, which is scaled as |A|j| oc J — J cl . Generally, an 
energy of the spin-triplet bound state is much lower than 
that of the spin-singlet bound stated 



C. Superconducting correlation 

For our model, we can consider four kinds of supercon- 
ducting correlations as mentioned in Sec lIIIl In order to 
estimate them numerically, we define the corresponding 
pair correlation functions as 



with 

Asi(r x ) = ^,1,1^,2,1 - ^,14^,24 (51) 

A T i(r x ) = ^,1,1^,2,1+^,1,4.^,2,1 (52) 

A S2 (7x) = ^,1,1^+1,24-^,14^+1,2,1 (53) 

A T2 (r x ) = ^,14^+1,24+^,14^+1,24, (54) 

which are calculated by the DMRG method. The cal- 
culated results for three sets of parameters are shown in 
Fig. For n = 0.825, U = 10, and J = 4 [Fig. 0[a)], 
the S2 pairing correlation is clearly the most dominant 
one, which shows a power-law length dependence. The 
ratio of the decay is estimated to be ~ r~ - 7 , which leads 
to = 0.714. This is consistent with the bosonization 
result for the spin gapful state. At n = 0.5, U = 10, 
J = 4 [Fig. 0b)], the system has a S = ground state 
but is somewhat closer to the fully-polarized ferromag- 
netic phase. The S2 pairing correlation is still the most 
dominant one but the Tl pairing correlation becomes 
much enhanced at a short distance \r x — r' x \ < 10. This 
is because the formation of a local spin-triplet bound 
state on a rung can gain some energies. We note that 
the T2 pairing correlation is also enhanced, reflecting a 
tendency to the fully-polarized ferromagnetic state. The 
decay ratio of the S2 correlation is ~ r -0 8 , which leads 
to A' + = 0.625. This value agrees well with our di- 
rect estimation of (See Sec lIV E[) . If the system fur- 
ther approaches the fully-polarized ferromagnetic phase 
[Fig. 0c)], the change in the correlation functions at 
short ranges becomes more prominent. Surprisingly, we 
find that the decay lengths of all the correlations are 
almost unchanged as far as J is fixed in the S = 
ground state. In the fully-polarized ferromagnetic phase 
[Fig. 0d)], the Tl pairing correlation is the most domi- 
nant one and only the S2 pairing correlation is competing. 
The decay ratio of both correlation functions is ~ 
which leads to K+ = 0.5. 



D. Anomalous Green function 

Let us now determine the pairing symmetry in the S2 
superconducting state. To this end, we calculate the one- 
particle anomalous Green function^ which exhibits the 
excitations of the Bogoliubov quasiparticles in the BCS 
theory. Therefore, we can see how the nodes appear in 
the superconducting gap function. The anomalous Green 
function is defined by 



Gk y ,k' y {Qx 



C 9x,fe H »T 



z-H + E 



V-'o 



(55) 



where 



^(K-^DHAj^A^)) 



ip™" ) denote the wavefunction of the ground state 



(50) with N e /2 up-spin and N e /2 down-spin electrons, and Eq 
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FIG. 9: (Color online) DMRG results for K p + as a function 
of J. The results for U = 10 and 20 at n — 0.5 are shown. 
Dotted line is guide for eyes. 



FIG. 

k„. = 



Left panel: 

N e /N 



Anomalous Green's function at k v 



0, 



5/8, U = 10, J = 2, and TV = 16. Right 



panel: Same as the left panel but at k y = ir, k' y — 0. 



relation Ffe 

-Fu..y (-k 



,k' u (k x ) = -Fk' y ,k y (k x ) and F ky ,k> y {k x ) 
). This indicates the formation of the d x 



wave like pairing state in our system, which is consistent 
with the weak-coupling results shown in Sec lIII Cl 



is chosen as the average value of E^(N^ — 1, — 1) and 
El (JVf , iVj ) . We then estimate the spectral function as 



F ky ,k' y (k x ) = --lmG ky ,k v (qx,u + ir)) (56) 
with rj = + and its frequency integral as 
F ky ,k y (k x ) = (^- 2 



(57) 



We calculate the anomalous Green function ([55)1 for the 
ladder with length L = 8 by using the Lanczos method. 
We should note that the BCS theory is not readily ap- 
plicable for (quasi) one-dimensional system; actually, 
Fk y ,k< (Qx) is not long-range ordered with a logarithmic 
decay as a function of L. However, apart from this pref- 
actor, we can naively expect that i*fe (q x ) calculated 
on finite-size systems should provide information on the 
pairing symmetry at intermediate distances. 

The calculated results for F ky ^(k x ,uj) are shown in 
Fig. [8j We can observe pronounced low-energy peaks 
at \k x \ = \ for all k y and k' y , which are the nearest 
to the Fermi momenta &f(= anc ^ much smaller 

peaks at higher energies for other momenta. We also 
find that the spectral weight vanishes at q x = 0, q x = tt, 
and k y — k'. The weights of the peaks appear to 
be similar to the BCS form of the condensation am- 
plitude, which has a maximum value at the Fermi mo- 
menta. We may thus assume that the superconduct- 
ing ground state in strongly correlated electron sys- 
tems can be characterized by a gap function, which is 
directly proportional to the frequency-integrated func- 
tion F kytk , y (q x )£2£L The function F ky . ky {k x ) obeys the 



E. Tomonaga-Luttinger liquid parameter 

For the calculation of the TL parameter, we use a re- 
cently proposed method of the DMRG technique.— As 
mentioned in Sec. IIII C[ the TL parameter deter- 
mines the long-range decays of the dry-like (S2) pair- 
ing and 4/cf-CDW correlations in the metallic TL-liquid 
ground state, whereas the parameter K~ is scaled as 
K~ — > 0. For the double-chain model, we define the 



TL parameter as 



K. 



2 g^O dq 



with 



N±(q) 



-Y, 



(58) 



(59) 



where 



n rx \ ± n rxt 2- In Fig. [51 we show the calcu- 



lated results for K p+ as a function of J for U = 10 and 20 
at n = 0.5. For [7=10, increases monotonously with 
increasing J, which is consistent with the monotonous 
increase in the binding energy with respect to J in the 
strong-coupling regime. For U — 20, the behaviors for 
J < 4 and J > 9 are quite similar to those for U = 10, 
although the values are somewhat smaller due to the ef- 
fects of the Umklapp scattering. However, we estimate 
the TL parameter as K+ = 0.5 for 5 < J < 8. This 
regime corresponds to the fully-polarized ferromagnetic 
phase and the TL parameter should be the same as that 
of a spinless fermion system. We thus find K p > 1/2 in 
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FIG. 10: (Color online) Phase diagram of our double-chain 
Hubbard model obtained from the DRMG calculations. The 
phase boundary between the fully-polarized ferromagnetic 
phase and d xy -\ike superconducting (SC) phase is shown for 
U = 5, 10, and 20. 

the entire region except in the fully-polarized ferromag- 
netic phase and we confirm that the S2-type supercon- 
ducting correlation is most dominant. 

F. Phase diagram 

In Fig. [TU1 we show the phase diagram of our model 
in the parameter space (n, J) for {7 = 5, 10, and 20. 
The phase boundary is determined by the calculated re- 
sults for the total-spin quantum number. At higher elec- 
tron concentrations from the phase boundary, the sys- 
tem is characterized as the S2-type spin-singlet super- 
conducting state, and at lower concentrations from the 
phase boundary, the system is characterized as the fully- 
polarized ferromagnetic state and simultaneously as the 
Tl-type spin-triplet superconducting state. We note that 
the fully-polarized ferromagnetic phase is spread to the 
higher concentration range with increasing U , and it oc- 
cupies the entire parameter region (n < 1, J > 0) at 
U = oo. This result is connected to the Nagaoka state- 
at U — ► oo. 



V. SUMMARY AND DISCUSSION 

We have studied the ground-state properties of the 
double-chain Hubbard model coupled with ferromagnetic 
exchange interaction by using the weak-coupling the- 
ory, DMRG technique, and the Lanczos method. We 



have thereby determined the ground-state phase diagram 
in the parameter space of the ferromagnetic exchange 
interaction and band filling. We have found that, in 
high electron density regime, the spin gap opens and 
the spin-singlet (i^y-wave-like pairing correlation is most 
dominant, and in low electron density regime, the fully- 
polarized ferromagnetic state is stabilized and simulta- 
neously the spin-triplet pj,-wave-like pairing correlation 
becomes most dominant. 

Here, let us make some comment on what happens 
if some additional terms are introduced to our model. 
First, we consider the case where the hopping integral 
between the two chains t± is added. In this case, the 
couplings g^jj and g£± become irrelevant, which leads 
to a collapse of the pairing state between different bands, 
i.e., k y = and ir. Thus, the spin-singlet superconduc- 
tivity is suppressed. Since the term t± induces the anti- 
ferromagnetic interaction on each rung, the spin-triplet 
superconductivity may also be suppressed. Next, we con- 
sider the case where the intersite Coulomb interaction 
between the two chains V± is added. In this case, we 
can easily imagine that the spin-triplet superconductivity 
is strongly suppressed because the cfy-wave-like pairing 
state is formed on rung. The double exchange interac- 
tion for the cizy-wave-like pairing state is also suppressed. 

Finally, let us discuss possible relationship between 
the present model and the model of two Hubbard chains 
coupled with zigzag bonds where the spin-triplet super- 
conductivity has been shown to occu n 17 ' 18 ! 19 The two 
models in the strong-coupling regime have the follow- 
ing features in common: (i) the superconductivity occurs 
near the region of ferromagnetism, (ii) there is a compe- 
tition between the spin-singlet and spin-triplet pairings, 
and (iii) the spin-gap is quite small or vanishes. These 
situations may well suggest that the spin-triplet pairing 
could be dominant in the present model near the phase 
boundary between the spin-singlet superconductivity and 
ferromagnetism. So far, we have not found any indica- 
tions that the spin-triplet pairing becomes more domi- 
nant than the spin-singlet pairing unless the ground state 
is spin polarized. We hope that this point will be clarified 
with more elaborate calculations done in future. 
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